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EXECUTIVE  SUMMARY 


This  report  was  sponsored  under  the  Federal  Aviation  Administration's  Aircraft 
Catastrophic  Failure  Prevention  Research  Program,  Mission  Need  Statement  No.  066- 
1 10.  This  report  documents  the  NASA  recommended  Matched-Filter-Based  One- 
Dimensional  Search  Method  for  the  gust  load  analysis  of  nonlinear  aircraft.  The  FAA 
Technical  Program  Monitor  was  Mr.  Terry  Barnes,  ANM-105N,  National  Resource 
Specialist,  Flight  Loads  and  Aeroelasticity;  FAA  COTR  was  Mr.  Thomas  DeFiore,  ACD- 
220,  Principal  Investigator,  Flight  Loads. 
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INTRODUCTION 


This  manual  and  the  accompanying  code  were  developed  as  partial  fulfillment  of  a 
NASA  agreement  with  the  FAA.  The  agreement  had  three  main  tasks.  First,  two 
NASA  developed  gust  analysis  methods  were  to  be  brought  to  the  same  level  of  maturity. 
These  analysis  methods  were  the  Matched-Filter-Based  (MFB)  method  and  the 
Stochastic-Simulation-Based  (SSB)  method.  Upon  completion  of  the  development  work, 
the  second  task  was  to  compare  the  methods  and  make  a  recommendation  selecting  which 
approach  was  best  suited  for  nonlinear  analyses.  This  work  was  completed  and  the 
results  given  in  a  presentation  at  the  Gust  Specialists  Meeting  in  LaJolla,  California  on 
April  22,  1993.  At  the  meeting,  it  was  recommended  by  NASA  that  the  MFB  one- 
dimensional  search  was  the  method  of  choice  for  time-correlated  gust  loads  analysis  of 
aircraft  with  nonlinear  systems.  The  third  and  final  task  was  to  develop  a  transportable 
computer  program  and  accompanying  documentation  for  using  the  recommended 
method. 

This  manual  describes  a  computer  code  called  "MFB  IDS"  which  performs  the  Matched - 
Filter-Based  one-dimensional  search.  The  MFB  one-dimensional  search  is  a 
deterministic  method  for  obtaining  maximized  and  time-correlated  design  loads  for 
aircraft  with  nonlinearities.  This  method  can,  however,  be  applied  to  both  linear  and 
nonlinear  aircraft.  The  paper  summarizes  the  method,  discusses  the  selection  of  gust 
intensity  for  the  method,  describes  the  FORTRAN  code  MFB  IDS  that  performs  the 
calculations,  and  presents  numerical  results  for  an  example  aircraft. 

BACKGROUND 

With  the  advent  of  aircraft  that  contain  large  numbers  of  nonlinearities  in  their  flight 
control  systems  and/or  gust  load  alleviation  systems,  existing  methods  for  certifying 
aircraft  for  gust  loads  may  not  be  adequate.  For  several  years  NASA  Langley  Research 
Center  has  conducted  research  in  the  area  of  time  correlated  gust  loads  and  has  published 
a  number  of  papers  on  the  subject  (refs.  1-6).  The  initial  research  was  restricted  to 
mathematically  linear  systems  (refs.  1-3).  Recently,  however,  the  focus  of  the  research 
has  been  on  defining  methods  that  will  compute  design  gust  loads  for  an  airplane  with  a 
nonlinear  control  system  (refs.  4-6).  To  date,  two  such  methods  have  been  defined:  one 
is  based  on  matched  filter  theory;  the  other  is  based  on  stochastic  simulation. 

The  Matched-Filter-Based  (MFB)  method  was  developed  first  and  was  reported  on  in 
reference  4.  The  MFB  method  employs  optimization  to  solve  for  its  answers  and  this 
method  comes  in  two  varieties;  the  first  uses  a  one-dimensional  search  procedure  and  the 
second  a  multi-dimensional  search  procedure.  The  first  is  significantly  fester  to  run  and, 
based  on  experience  with  a  number  of  nonlinear  models  gives  design  loads  only  slightly 
lower  in  magnitude  than  the  second. 
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The  Stochastic-Simulation-Based  (SSB)  method  has  evolved  over  the  past  two  years. 
References  5  and  6  describe  the  method.  In  reference  6  a  comparison  of  the  MFB  and 
SSB  methods  was  made.  The  results  predicted  by  the  two  methods  are  strikingly  similar 
and  demonstrate  that  the  key  quantities  from  the  MFB  method  (viz.  critical  gust  profile, 
maximized  load,  and  time-correlated  load)  are  realizable  in  a  stochastic  analysis. 

Another  significant  finding  in  inference  6  was  the  relative  computational  costs  of 
performing  analyses  using  ti  e  MFB  and  SSB  methods.  For  linear  models,  the  MFB 
method  is  much  more  efficient  than  the  SSB  method.  For  nonlinear  models,  the  MFB 
multi-dimensional  search  is  much  more  expensive  that  the  SSB  method,  while  the  MFB 
one-dimensional  search  requires  less  time  that  the  SSB  method. 

Since  the  MFB  multi-dimensional  search  method  as  now  implemented  is  prohibitively 
expensive,  the  options  for  methods  that  can  practically  he  applied  to  nonlinear  systems 
are  the  MFB  one-dimensional  search  and  the  SSB  methods. 

Based  on  the  results  in  reference  6,  the  one-dimensional  search  is  able  to  predict  the 
maximized  loads  for  nonlinear  systems.  In  addition,  it  requires  less  computer  resources 
than  the  SSB  method.  These  factors  show  the  utility  of  u.sing  the  MFB  one-dimensional 
search  as  a  means  of  obtaining  time-correlated  gust  loads  for  aircraft  with  nonlinear 
control  systems. 

REVIEW  OF  THE  MFB  METHOD  FOR  LINEAR  SYSTEMS 

The  purpose  of  this  section  of  the  manual  is  to  review  the  basic  matched  filter  concepts. 
A  detailed  theoretical  development  of  the  MFB  linear  method  can  be  found  in  reference 
2.  The  signal  flow  diagram  in  figure  1  outlines  the  implementation  and  illustrates  the 
intermediate  and  final  products  of  the  method. 

Transfer-function  representations  of  atmospheric  turbulence  and  airplane  loads  are 
combined  in  series  and  repre.sent  the  "known  dynamics"  boxes  in  the  figure.  A  transfer- 
function  representation  of  the  von  Karman  spectrum  in  chosen  for  the  gust  filter.  Load  y 
is  the  load  to  be  maximized.  Loads  z  j  through  Zn  are  the  loads  to  be  time  correlated  with 
load  y.  There  are  three  major  steps  in  the  process: 

Step  i  The  application  of  an  impulse  function  of  unit  strength  to  the 

combined  linear  system,  producing  the  impulse  response  of  load  y, 
h(t). 

Step  ii  The  normalization  of  this  impulse  respon.se  by  the  square  root  of  its 
energy  and  then  reversing  it  in  time. 

Step  iii  The  application  of  this  normalized  reversed  signal  to  the  combined 
linear  system,  producing  time  histories  of  load  y  and  time  histories 
of  loads  zi  through  Zi,. 

For  simplicity  of  discussion  these  three  steps  will  be  referred  to  as  the  "MFB  linear 
method." 
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The  square  root  of  the  energy  is  defined  by 


where  h(t)  is  the  impulse  response  of  load  y.  For  a  stable  system  the  numerical  value  of 
the  quantity  on  the  right  side  of  equation  (1)  approaches  a  constant  value  as  to  is 
increased. 

The  selection  of  to  is  based  on  the  time  required  for  the  load  impulse  responses  to  damp 
out  to  near  zero.  Figure  1  shows  that  the  impulse  response  of  load  y  has  essentially 
damped  out  at  time  to-  Thus,  to  is  large  enough  for  the  response  shown  in  figure  1 .  Too 
large  a  to  value  will,  however,  unduly  increase  the  amount  of  computations  required.  The 
analyst  must  choose  a  value  that  provides  accurate  results  while  minimizing  the  use  of 
computer  resources. 

Within  the  time  history  of  load  y  in  step  iii,  the  maximum  value  is  ymax-  For  linear 
systems,  the  theory  guarantees  that  no  other  signal  similarly  normalized  will  produce  a 
value  of  y  larger  than  ymax-  This  guarantee  is  a  fundamental  result  of  the  Mra  linear 
method. 

MFB  ONE-DIMENSIONAL  SEARCH  DESCRIPTION 

The  goal  of  Matched  Filter  Theory  as  applied  to  nonlinear  systems  is  the  same  as  that  for 
linear  systems:  to  find  the  maximized  response  time  history,  the  maximum  value  of  the 
response  within  that  time  history,  and  the  time-correlated  response  time  histories. 

Because  the  systems  are  not  linear,  the  superposition  principle  of  the  MFB  linear  method 
no  longer  holds  and  the  solutions  for  maximized  loads  cannot  conveniently  be  obtained 
directly.  The  only  practical  means  of  finding  the  excitation  waveform  that  maximizes 
ymax  is  a  search  procedure.  The  search  is  conducted  systematically,  subject  to  the 
constraint  that  the  excitation  wavefoiTn  have  a  "unit"  energy. 

Because  superposition  no  longer  holds,  the  magnitude  and  character  of  the  responses  are 
not  necessarily  proportional  to  the  magnitude  of  the  input.  For  the  remainder  of  this 
paper  two  input  magnitudes  are  important:  k,  the  strength  of  the  initial  impulse;  and  Og, 
the  design  value  of  the  gust  intensity.  (For  the  MFB  linear  method,  the  magnitude  of 
both  of  these  quantities  was  unity.)  For  nonlinear  systems  and  a  specific  og,  the  shape  of 
the  excitation  waveform  is  a  function  of  k  ,  and,  consequently,  the  quantity  y^ax  is  also  a 
function  of  this  parameter. 

The  one-dimensional  search  procedure  performs  a  systematic  variation  of  the  quantity  k 
to  find  the  shape  of  an  excitation  waveform  that  maximizes  ymax-  Figure  2  contains  a 
signal  flow  diagram  for  this  search  procedure.  Figure  2  is  very  similar  to  figure  1,  but 
contains  some  subtle  yet  important  differences  that  are  indicated  by  the  shaded  boxes  and 
by  quotation  marks.  In  figure  2  the  initial  impulse  has  a  non-unity  strength;  the  aircraft 
loads  portion  of  the  known  dynamics  box  contains  nonlinearities;  and  the  shape  of  the 
excitation  waveform  and  the  value  of  ymax  are  functions  of  the  initial  impulse  strength. 

In  addition,  the  "matched"  excitation  waveform  and  the  "matched"  load  are  shown  in 
quotes  because,  for  nonlinear  .systems,  there  is  no  guarantee  that  ymax  is  a  global 
maximum. 
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The  application  of  the  one-dimensional  search  procedure  is  as  follows: 


Step  1  Select  a  specific  design  value  for  ag. 

Step  2  Select  a  set  of  values  for  k. 

Step  3  For  each  value  of  k,  perform  steps  i  through  iii  of  the  MFB  linear 
method  to  obtain  a  set  of  "ymax"  values  and  also  the  corresponding 
"matched"  excitation  waveforms. 

Step  4  From  step  3  above,  find  the  maximum  value  of  ymax  ^tnd  its 
corresponding  "matched"  excitation  waveform. 

As  seen  in  the  first  pass  (upper  half)  of  figure  2,  the  variation  of  the  impulse  strength  (k) 
affects  the  MFB  one-dimensional  search  analysis  by  changing  the  shape  of  the  excitation 
waveform.  For  sufficiently  low  impulse  strengths,  the  shape  of  the  excitation  waveform 
nonlinear  models  will  be  invariant  with  k.  While  in  this  invariant  region,  the 
excitation  waveforms  will  be  similar  to  those  that  are  obtained  for  linear  models.  For 
larger  intensities  the  system  nonlinearities  are  engaged  and  will  cause  the  impulse 
responses  and  corresponding  excitation  waveforms  to  change  shape.  Consequently,  they 
will  differ  from  those  obtained  by  linear  systems. 

As  seen  in  the  second  pass  (lower  halO  of  figure  2,  the  gust  intensity  affects  the  MFB 
one-dimensional  search  analysis  by  scaling  the  excitation  waveform  prior  to  being 
applied  to  the  nonlinear  model.  Consequently,  a  low  gust  intensity  should  result  in  the 
nonlinear  model  behaving  linearly.  As  gust  intensity  is  increased  beyond  some  threshold 
the  nonlinear  model  response  will  begin  to  deviate  from  that  of  its  linear  counterpart. 

The  effects  of  varying  impulse  strength  and  gust  intensity  will  be  discussed  further  in  the 
numerical  results  section  of  the  paper. 

SELECTION  OF  GUST  INTENSITIES 

The  purpose  of  this  section  of  the  paper  is  to  present  the  reasoning  behind  the  selection  of 
the  values  of  og.  Reference  6  describes  the  selection  of  gust  intensities  for  the  MFB  and 
SSB  methods. 

The  following  equation,  from  reference  7,  expresses  the  "design  value"  of  quantity  y  as 
defined  in  the  design  envelope  criterion 

(2) 

where  the  quantity  Ay  is  the  RMS  value  of  quantity  y  per  unit  RMS  gust  intensity, 
obtained  from  a  conventional  random  process  analysis  of  the  airplane  and  Uo  is  specified 
in  the  criterion.  From  reference  8  the  quantity  Uo  in  equation  (2)  is  shown  to  be  the 
product  of  the  gust  RMS  value  and  the  design  ratio  of  peak  value  of  load  to  RMS  value  of 
load,  and,  therefore  the  quantity  yaesign  is  interpreted  as  a  peak  value. 
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Reference  4  shows  that,  as  a  consequence  of  the  normalization  of  the  excitation 
waveform  by  the  square  root  or  its  own  energy  and  the  use  of  unity  gust  intensity,  the 
quantity  ymax  from  the  MFB  linear  method  is  equal  to  the  quantity  a,  from  a 
conventional  random  process  analysis,  or 

In  equation  (3)  ymax  is  interpreted  as  an  RMS  value,  not  a  peak  value.  Substituting 
equation  (3)  into  equation  (2),  yjesign  is  now 

y design  ~  ymax(^g  ~  ^ 

If,  in  performing  the  MFB  linear  method,  Uo  is  used  for  the  gust  intensity  then  the 
quantity  ymax  is  equal  to 

ymax^^g  ~  ~ 

The  right  hand  sides  of  equations  (2)  and  (5)  are  seen  to  be  equal,  therefore 

y design  ~  ymax^^g  ~ 

Two  options  for  the  value  of  og  have  been  offered:  Og  =  1,  for  which  ydesign  is  defined  by 
equation  (4);  and  og  =  Uo,  for  which  yue.sign  is  definea  by  equation  (6).  When  analyzing  a 
linear  system  the  choice  of  og  is  inelevant  because  the  same  value  of  ydesign  will  be 
obtained  in  either  case.  However,  when  nonlinearities  are  introduced  into  aircraft  control 
systems,  loads  are  not  simply  proportional  to  gust  intensity.  Consequently,  Og  should  be 
set  to  Uo  in  the  MFB  nonlinear  calculations,  or 

^sMFB  ~  ^ a 

and  the  resulting  "ymax"  values  from  the  method  should  be  interpreted  as  ydesign- 


DESCRIPTION  OF  MFBIDS 

MFB  IDS  is  a  FORTRAN  77  program  which  performs  the  Matched-Filter-Based  one¬ 
dimensional  search.  The  MFB  IDS  solution  procedure  follows  what  was  outlined  in 
figure  2.  The  code  must  be  run  once  for  each  combination  of  gust  intensity  and  output 
quantity  for  which  a  maximized  value  is  desired. 

Figure  3  shows  the  solution  procedure  used  by  MFBIDS.  The  program  first  reads  the 
input  data  then  generates  the  impul.se  responses  for  each  k  value  by  calling  the  simulation 
subroutine.  The  simulation  subroutine  uses  a  public  domain  ordinary  differential 
equation  solver  to  generate  output  time  histories  using  a  user  defined  subroutine 
containing  the  aircraft  equations  of  motion.  Next,  the  excitation  waveforms  are  generated 
by  reversing  the  impulse  responses  in  time  and  normalizing  them  by  the  square  root  their 
respective  energies.  Then,  the  simulation  subroutine  is  again  called  for  each  of  the 
normalized  excitation  wavefonns  to  obtain  the  "maximized"  load  response  time  histories. 
Finally,  the  pertinent  output  quantities  are  written  to  computer  files. 


This  section  of  the  manual  describes  the  main  parts  of  MFB  IDS. 

Required  Files 

Six  files  are  required  to  run  MFB  IDS. 

MFBIDS.F  Main  program,  see  Appendix  A 

MFB  IDS. INC  File  containing  the  common  blocks  used  by 

MFB  IDS. F,  see  Appendix  B 

MODEL. F  User  .supplied  file  containing  the 

subroutine  EQSMOT,  the  aircraft  equations  of 
motion,  see  Appendix  C 

MODEL.INP  User  supplied  file  containing  the  input 

quantities  required  to  run  the  code,  see  Appendix  D 

LSODE.F  Public  domain  ordinary  differential  equation 

solver 

INTUTILS.F  Public  domain  subroutines  used  by  LSODE.F 

The  two  public  domain  files  LSODE.F  and  INTUTILS.F  are  well  documented  in  their 
respective  source  codes  and  will  not  be  discussed  in  detail  in  this  manual.  These  files 
were  selected  because  of  the  fact  that  they  are  in  the  public  domain.  The  user  is 
encouraged  to  substitute  more  efficient  ordinary  differential  equation  solvers  if  available. 

MFBIDS.F  Subroutines 

Subroutine  DATAIN  Reads  the  input  parameters  from  the  input  file 

Subroutine  SIMULATE  Computes  output  time  histories  using  LSODE.F 

and  .MODEL.F 

Subroutine  SAVMATFORM  Saves  output  quantities  to  a  MATRIX;  ^  readable 

file 

Subroutine  MATSAV  Used  by  Subroutine  SAVMATFORM  to  save  data 

in  MATRIX  X  format 

Subroutine  SAVASCIIFORM  Saves  output  quantities  in  a  more  easily  readable 

ascii  file 

MODEL.F 

MODEL.F  is  a  user  supplied  file  containing  the  subroutine  EQSMOT,  the  equations  for 
the  gust  filter  in  series  with  the  aircraft  equations  of  motion.  This  subroutine  uses 
variables  x  (vector  of  states)  and  u  (excitation)  to  obtain  xd  (derivative  of  x)  and  y  (vector 
of  outputs). 
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The  following  lines  of  code  must  be  present  at  the  top  of  the  file  MODEL.F: 

subroutine  EQSMOT(neq,t,x.xd) 
parameter(maxout=2()) 
double  precision  x(*),  xd(*),  y(maxout) 
double  precision  t,  tstart,  tend,  uO,  ul 
common  /eqsmotcom/y.tstart.tend.uO.u  1 
c 

u=(t-tstart)*(ul-uO)/(tend-tstart)+uO 

In  spite  of  the  fact  that  the  user  specifies  a  desired  time  step,  the  ordinary  differential 
equation  solver  LSODE  calls  the  EQSMOT  subroutine  at  values  of  time  (t)  between  time 
steps.  Consequently,  the  equation  on  the  last  line  of  the  above  code  is  included  to 
provide  a  linear  interpolation  for  values  of  u  between  time  steps. 

A  linear  system  can  be  u.sed  to  demonstrate  how  x,  xd,  y  and  u  are  related.  For  a  linear 
system  only,  the  equations  of  motion  can  be  written  in  state  space  form; 

{xd}  =  [A]{x}  +  [B]u 
{y}  =  [C]{x}  +  [D]u 

While  equation  (8)  is  linear,  the  relationship  between  x,  u  and  xd  is,  in  general,  a 
nonlinear  one. 

Description  Of  Impulse  Input 

Since  the  impulse  generating  procedure  is  one  of  the  key  components  of  the  program  and 
there  are  several  methods  that  can  be  used  to  generate  impulse  functions,  a  few  words 
describing  the  impulse  function  procedure  implemented  in  the  program  are  warranted. 
The  procedure  chosen  to  generate  the  impulse  function  in  this  code  is  straightforeword 
and  can  be  found  in  main  program  listing  in  Appendix  A.  In  general,  the  impulse  input 
function,  as  seen  by  the  differential  equation  solver,  is  a  ramp  up  from  zero  at  the  first 
time  step  to  a  value  of  k/(2*deltat)  at  the  second  time  step,  a  constant  value  of 
k/(2*deltat)  between  time  steps  two  and  three,  and  a  ramp  down  to  zero  between  time 
steps  three  and  four. 

•  Description  Of  Output 

Key  information  is  written  to  standard  out  (unit  6)  as  shown  in  Appendix  E.  This 

*  information  allows  the  user  to  monitor  the  progress  of  the  program.  Two  additional 
output  files  can  be  created  which  contain  time  history  and  correlated  output  information. 
Both  of  these  output  files  are  ascii  files:  one  is  in  a  MATRIXx  readable  form  and  the 
other  is  in  a  more  easily  read  form.  The  variables  contained  in  the  MATRIXx  readable 
file  are  discussed  below.  The  quantities  found  in  the  other  output  file  are  self  explanatory 
and  will  not  be  discus.sed  here. 

The  following  scalar  quantities  are  written  to  the  MATRIXx  file: 
sigmag  Gust  intensity, 

noutmx  The  output  quantity  to  be  maximized. 
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deltat 


The  lime  step. 

tmaximp  The  length  of  the  impulse  respon.ses. 

The  following  arrays  can  also  written  to  the  MATRIX x  file: 

allkvals  This  vector  has  length  nkvals,  and  stores  all  the  k  values  used  in 
the  analysis. 

maxout  The  maximum  output  values  for  each  k  value  are  stored  in  each 

column.  Thus,  this  array  will  be  an  nout  by  nkvals  array. 

impres#  An  anay  of  this  form  is  created  for  each  output  quantity.  For 

in.stance,  for  output  1  an  array  named  "impres  1"  will  be  created  and 
for  output  5  an  airay  named  "impresS"  will  be  created.  Each 
column  of  these  arrays  are  impulse  response  time  histories  for  one 
of  the  k  values.  Thus,  these  are  ntsteps  by  nkvalues  arrays. 

wavef#  An  anay  of  this  form  is  created  for  the  output  to  be  maximized. 

For  instance,  if  output  6  is  to  be  maximized,  then  an  array  named 
"wavef6"  will  be  created.  Each  column  of  this  array  is  an 
excitation  waveform  matched  to  output  #  for  a  k  value.  Thus,  this 
is  an  ntsteps  by  nkvalues  array. 

exresp#  An  anay  of  this  form  is  created  for  each  output  quantity.  For 
instance,  for  output  1  an  array  named  "exresp  1"  will  be  created. 
Each  column  of  this  array  is  an  excitation  waveform  response  for  a 
k  value.  Thus,  these  are  (2*ntsteps+l)  by  nkvalues  arrays. 

NUMERICAL  EXAMPLE 

This  section  describes  the  step  by  step  process  of  obtaining  numerical  results  at  a  gust 
intensity  of  1,530  in/sec. 

Step  1  -  Create  a  subroutine  containing  the  gust  filter  in  series  with  the  aircraft 
equations  of  motion. 

The  nonlinear  simulation  model  of  the  ARW-2  drone  aircraft  equipped  with  a  nonlinear 
control  system  was  used  in  this  example.  This  model  was  connected  in  series  with  a 
transfer-function  representation  of  atmo.spheric  turbulence.  The  transfer  function  used  in 
this  example  was  (ref.  8) 

pT  [1  +  2.6]8(L/V)s][l  -h  0.1298(L/V)s] 

T]  [1  +  2. 083(LIV)s ][I  +  0. 823(L/V)s ][1  +  0. 0898(LIV)s ] 

which  approximates  the  square  root  of  the  von  Karmon  power  spectral  density  function. 
The  quantity  Og  is  the  intensity  of  the  gust  or  standard  deviation  —  which,  assuming  zero 
mean,  is  also  equal  to  the  root-mean-.square,  or  RMS,  value  —  of  gust  velocity. 

Figure  4  shows  a  block  diagram  of  the  ARW-2  simulation  model  and  includes  the 
aeroelastic  plant,  a  gust  load  alleviation  (GLA)  control  law,  and  nonlinear  control 
elements.  The  aeroelastic  plant  is  a  linear,  s-plane  aeroelastic  half-model  consisting  of 
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two  longitudinal  rigid-body  modes  and  three  symmetric  flexible  modes.  Unsteady 
aerodynamics  were  obtained  using  the  doublet  lattice  method  (ref.  9).  The  model  also 
includes  the  dynamics  of  the  control  surface  actuators.  The  two-input/two-output  GLA 
control  law  was  obtained  using  a  Linear  Quadratic  Gaussian  design  approach  with  the 
intent  of  reducing  wing  root  bending  moment  (ref.  10).  The  nonlinear  elements  impose 
deflection  limits  of  ±1”  on  the  elevators  and  0“  to  +1“  on  the  ailerons  to  simulate  spoilers 
The  model  contains  32  states,  and  the  analysis  conditions  are  at  a  Mach  number  of  0.86 
and  an  altitude  of  24, (XK)  feet. 

The  resulting  combined  system  has  a  single  input  (white  noise)  and  many  outputs 
including  the  gust  velocity  (y(14)).  The  output  quantities  for  this  model  are  as  follows: 


y(i)-y(3) 

=  No  physical  interpretation 

y(4) 

=  Tip  acceleration 

y(5) 

=  Fuselage  acceleration 

y(6) 

=  Wing  root  bending  moment  (WRBM) 

y(7) 

=  Wing  root  shear  (WRS) 

y(8) 

=  Wing  outboard  bending  moment  (WOBBM) 

y(9) 

=  Wing  outboard  torsion  moment  (WOBTM) 

y(i0) 

=  Elevator  deflection 

y(ii) 

=  Elevator  rate 

y(i2) 

=  Aileron  deflection 

y(i3) 

=  Aileron  rate 

y(i4) 

=  Gust  velocity 

Appendix  C  contains  the  FORTRAN  version  of  the  analytical  model.  This  model  was 
created  using  MATRIXx  SYSTEM  BUILD  (ref.  1 1)  and  converted  to  FORTRAN  using 
the  HYPERCODE  (ref.  12)  package.  As  shown  in  Appendix  C,  the  appropriate  lines  of 
code  were  added  to  the  model  created  by  the  HYPERCODE  package  to  make  it 
compatible  with  the  differential  equation  solving  scheme  used  in  MFBIDS. 

Step  2  -  Create  an  input  file  for  the  model. 

The  length  of  the  state  vector  and  number  of  model  output  quantities  are  contained  in  the 
input  file  MODEL.INP  and  must  be  made  compatible  with  the  analytical  model  created 
in  step  1.  Other  input  quantities  like  the  gust  intensity,  identifying  the  output  quantity  to 
be  maximized,  and  the  range  of  k  values  must  also  be  selected.  (Trial  and  error  will  be 
required  to  find  the  range  of  k  values  where  the  maximum  load  value  is  obtained.)  In 
addition  to  these  parameters,  the  length  of  time  for  the  impulse  responses  and  the  time 
step  must  be  chosen.  The  length  of  the  impulse  responses  should  be  only  long  enough  to 
allow  the  output  quantity  to  damp  out  to  a  relatively  small  value  as  shown  in  the  example 
in  figure  4  (a),  while  the  time  step  should  be  chosen  in  accordance  with  the  model 
response  characteristics. 

The  input  file  used  in  the  ARW-2  analysis  is  shown  in  Appendix  D. 

Step  3  -  Modify  the  parameters. 

Appendix  B  contains  a  listing  of  MFB  IDS.INC  and  a  description  of  the  parameters  and 
variables  contained  there.  The  parameters  found  in  the  file  MFB  IDS.INC  and  MODEL.F 
must  be  made  large  enough  to  accommodate  the  analytical  model  (step  1)  and  the  input 
file  (step  2).  Computer  storage  can  be  minimized  by  making  parameters  "maxstates"  and 
"maxout"  equal  to  the  minimum  values  required  for  a  given  analytical  model. 
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Sicp  4  -  Compile  llie  code. 

The  following  stateincnl  compiles  Jhe  cotie  on  a  SUN  workstation. 
m  MPB  IDS.P  MODEL.F  LSODE.F  INTUTfLS.F 


Slop  5  -  Execiifc  llie  code. 

Appcndi,\  D  shows  the  initial  input  file  used  in  the  analysis  of  the  ARW-2  model.  The 
initial  value  of  gust  intensity  was  l.5.^0  in/sec  (1.5X85  ft/scc).  This  initial  value  of  gust 
intensity  was  chosen  so  that  the  nonlineaiities  would  be  inv(>ked  in  this  example.  Nine  k 
values  were  used  ranging  from  10  to  I5,(XK).  The  output  cpiantity  to  be  maximized  was 
y(b).  wing  root  bending  moment.  VVRBM.  The  standard  output  from  this  run  is  shown  in 
Appendix  E. 

The  following  statement  w  ill  run  the  code  on  a  UNIX  system 
a.out  <  MODEL.INP 

Step  6  -  Kxnniinc  output. 

Of  the  nine  k  values  used  in  the  analysis,  all  the  intermediate  results  for  three  of  the  nine 
k  values  are  shown  iti  figure  5.  Beginning  with  plot  (a),  the  impul.se  respon.ses  for  the 
three  k  values  arc  shown.  Plot  (b)  gives  the  corresponding  excitation  waveforms  obtained 
from  these  impulse  responses.  Plot  (c)  .shows  the  critical  gust  profiles  and  plot  (d)  shows 
the  loaximizcd  load  responses  created  by  the  gust  input.  Plot  (e)  depicts  maximum 
VVRBM  as  a  function  of  initial  impulse  strength:  the  circles  are  actual  numerical  results; 
the  solid  line  is  a  faired  line.  For  this  particular  example,  a  value  of  k  of  about  2.410 
creates  an  excitation  wavcF'ini  and  critical  giisT  profile  that  yields  the  largest  maximized 
value  of  VVRBM.  296.gy  4  in-lbs. 

Step  7  -  Vary  |)aramcter,s. 

Flic  user  may  wish  to  refine  the  results  shown  in  figure  5.  With  this  in  mind,  a  new  range 
of  k  values  w  .-is  selected  and  the  case  was  rerun.  Figure  6  shows  the  results  of  using  the 
new  set  <if  k  values  distributed  between  4(K)  and  6.000  producing  a  much  smoother  curve. 
The  maximum  load  value  obtained  using  the  refined  k  distribution  was  296,804  in-lbs  at  a 
k  value  of  approximately  2,1  ly.  For  this  particular  airplane  and  set  of  flight  conditions  a 
larger  load  \  aluc  was  not  obtained  with  the  new  range  of  k  values. 

EXAMINATION  OF  RFiJUMS 

This  section  has  been  incliulcd  in  the  manual  to  provide  the  user  with  insight  into  how  the 
results  Irom  Ml  B IDS  should  be  interpreted.  Ihc  results  presented  in  the  previous 
.section  will  be  shown  along  \j[ith  results  obtained  using  two  additional  gust  intensities. 

To  obtain  results  at  different  gust  intensities,  the  sigma  value  of  the  input  file  shown  in 
Appendix  D  was  changed  and  MFBIDS  wa.s  rerun.  In  addition,  answers  were  obtained 
for  each  of  Ihc  three  gust  intensities  using  a  linearized  version  of  the  model.  For 
comparison  these  linear  answci’s  will  be  plotted  with  the  results  from  Ihc  nonlinear 
model.  Note  Ih.'it  when  using  a  linear  model,  only  one  k  \  alue  needs  to  be  used  because 
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the  answer  will  not  be  a  function  of  k.  Consequently,  the  maximum  load  value  plotted 
versus  k  is  a  horizontal  line  for  a  linear  model. 

The  results  for  gust  intensity  values  of  1,020  in/sec,  1,530  in/sec  and  2,040  in/sec  are 
shown  in  figure  7.  Before  proceeding,  it  is  important  to  point  out  that  results  are  problem 
dependent,  and  while  these  results  are  typical  of  all  the  aircraft  examined  by  NASA, 
different  aircraft  with  different  nonlinearities  may  exhibit  different  trends.  With  this  in 
mind,  two  important  points  will  be  made  concerning  the  general  trends  exhibited  in  figure 
7.  First,  regardless  of  the  gust  intensity  value,  the  maximum  load  is  constant  for  small 
values  of  k.  Consequently,  when  searching  for  the  maximum  load,  the  range  of  interest 
for  the  k  values  can  be  limited  on  the  low  end.  Second,  the  maximum  load  value 
decreases  toward  some  relatively  small  value  as  k  is  made  very  large.  Thus,  the  range  of 
interest  for  k  values  can  also  be  limited  on  the  high  end.  What  happens  in  between  these 
two  extremes  depends  on  the  specific  gust  intensity. 

In  figure  7  the  maximum  load  values  are  constant  for  k  less  that  400.  In  plot  7(a)  a 
significantly  larger  load  value  was  not  found  when  k  was  increased  beyond  400,  while  for 
plots  7(b)  and  7(c)  a  larger  value  was  found.  This  indicates  that  the  character  of  the 
results  (i.e.,  the  shape  of  the  maximum  load  versus  k  plot)  and  the  specific  k  value  that 
produces  the  maximum  load  value  are  functions  of  gust  intensity.  But,  the  range  of  k 
values  where  the  maximum  will  be  found  is  generally  limited,  and  in  this  case  the  range 
is  roughly  between  400  and  10,000.  It  should  also  be  noted  that  the  difference  between 
answers  from  linear  and  nonlinear  models  is  also  a  function  of  the  gust  intensity.  Larger 
gust  intensity  values  generally  result  in  larger  differences  between  answers  from  linear 
and  nonlinear  models.  Here,  this  difference  is  2%  for  the  lowest  gust  intensity  and  18% 
for  the  largest  gust  intensity. 


CONCLUDING  REMARKS 

This  manual  has  reviewed  the  theory  behind  the  Matched-Filter-Based  one-dimensional 
search  procedure.  The  code  that  performs  this  procedure  has  been  discussed  and  example 
numerical  results  were  presented  and  interpreted.  The  code,  MFBIDS,  is  available  in  a 
self  contained  form.  It  has  all  the  required  equation  solvers  and  files  necessary  to  run  the 
example  problem.  The  user  is,  however,  encouraged  to  modify  the  existing  code  by 
inserting  more  efficient  routines  if  available. 
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Figure  1 .  MFB  linear  method  signal  flow  diagram. 
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Figure  2.  Nonlinear  MFB  signal  flow  diagram  for  the  one-dimensional  search. 
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Figure  3.  MFB IDS  solution  procedure. 
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Figure  4.  Nonlinear  block  diagram  for  the  ARW-2  drone  aircraft. 
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(a)  Og=  1,020  in/sec  (85  ft/sec) 


ICO  '000  10000 


k 

(b)  (Tg  =  1 .530  In/sec  (1 .5X85  ft/sec) 


<1  too  1000  10000 

k 


f c)  Og  =  2,040  in/sec  (2X85  ft/sec) 

Figure  7.-  Maximum  WRBM  versus  k  for  gust  intensities. 
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Appendix  A  •  MFBIDS.F  source  code  listing 


program  MFBIDS 

,*i(c  +  +  +  +  +  *  +  +  **  +  +  **  +  *  +  *  +  +  +  +  +  +  >(‘*******<l‘******  +  *****  +  +>f‘  +  +  +  +  >(<***  +  * 


c*  MFBIDS.F; 


c*  written  by; 


a  matched- filter-based  method  of  obtaining 
maximized  and  time-coirelated  gust  loads  for 
linear  and  nonlinear  aircraft 

Robert  C.  Scott 

NASA  Langley  Research  Center  MS/34() 
Hampton,  VA  23665 

804-864-2838 


* 

* 

* 

* 

* 

* 

* 

* 

* 

* 


c*  compile:  til  MFB 1  DS.F  MODEL.F  LSODE.F  INTUTILS.F  * 

c*  * 

c*  where  MODEL.F  is  the  a/c  model  * 


include  'MFBIDS.INC 

character*  80  ctitle 

character*  80  cease 

common  /case/ccase 

common  /title/ctitle 

dimension  yth(maxout,2*maxtsteps) 


c  **read  input  parameters 
call  DATAIN 


write(6,*)  ctitle 
write(6,*)  cease 

c  **calculate  number  of  time  steps  for  impulse  response 
ntsteps=int(tmaximp/deltat  -t-  ().(M)1)  +  1 
write(6,*) '  ’ 

write(6,*) '  length  of  simulation  for  impulse  responses  = 

+  tmaximp 

write(6,*) '  number  of  time  steps  for  impulse  responses  = 

•f  ntsteps 
write(6,*) ' ' 

write(6,*) '  length  of  simulation  for  excitation  responses  = 
+  2*tmaximp 

write(6,*) '  number  of  time  steps  for  excitation  responses  = 
+  2*ntsteps-l 

write(6,*) ' ' 

write(6,*) '  gust  intensity  =  sigma 
write(6,*) ' ' 

c  **create  array  of  zeros  for  inpulse  input 
do  100  i=l,(2*ntsteps-l) 

100  aimpt(i)^.0 

c  **calculate  k  values 
do  200  k=Lnkvals 
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if  (nkvals.eq.l)  then 
allkvals(l)=akmin 
else 

allkvals(k)=  10.0’'‘*(  (k- 1  )*(log  10(akinax)-log  lO(akmin))/ 

+  float(nk vals- 1  )+log  1  ()(akmin) ) 

endif 

200  continue 

c  **obtain  impulse  response  for  each  k  value 

write(6,*)  'calculating  impulse  responses  for  each  k  value' 
do  300  k=l,nkvals 
ak=allkvals(k) 

c  **assign  impulse  strength  k  to  array  aimpt 
aimpt(2)=ak/deltat/2 
aimpt(3)=ak/deltat/2 

c  **call  subroutine  SIMULATE  to  obtain  impulse  response 
call  SIMULATE(yth,aimpt,ntsteps) 
do  300  iout=l,nout 
do  300  j=l,ntsteps 

300  aimp_res(j,iout,k)=yth(iout,j) 

c  **calculate  normalized  excitation  waveforms 

write(6,*)  'calculating  normalized  excitation  wavefonns' 
do  400  k=l,nkvals 
i=noutmx 
energy=0.0 

energy=aimp_res(l,i,k)**2  +  aimp_res(ntsteps,i,k)**2 
do  401  j=2,ntsteps-i 

401  energy  =  energy+2*(aimp_res(j,i,k))**2 

energy  =  sqrt(  tmaximp*energy*.5/float(ntsteps)/3. 14159  ) 
write(6,*) '  ',  k, '  k  =  ',  allkvals(k), 

+  '  sqrt(energy)  of  output',  noutmx, '  = ',  energy 

do  402  j=ntsteps.  1 1 

402  wave(ntsteps- j+ 1  ,k)=sigma*aimp_res(j,i,k)/energy 
400  continue 

c  **obtain  maximized  responses 

c  Note  that  here  the  length  of  the  simulation  is  2*tmaximp  and  the 
c  number  of  time  steps  for  maximized  responses  is  2*ntsteps- 1 

write(6,*)  'calculating  maximized  responses' 
do  500  k=l,nkvals 

do  501  i=l,ntsteps 

501  aimpt(j)=wave(j,k) 

do  502  j=(ntsteps+l),  (2*ntsteps-l) 

502  aimpt(j)=0.0 

call  SIMULATE(yth,aimpt,  (2*nts*eps-l) ) 
do  503  j=l,(2*ntsteps-l) 
do  503  i=I,nout 

503  aexc_res(j,i,k)=yth(i,j) 


write(6,*) '  k,'  k  =  ',allkvals(k), 

+  '  maximum  value  of  output 

+  noutmx, '  = aexc_res(ntsteps,noutmx,k) 

500  continue 

c  **save  data 

write(6,*)  'saving  data’ 

if  (iouttype.eq.l.or.iouttype.eq.3)  call  SAVASCIIFORM 
if  (iouttype.eq.2.or.iouttype.eq.3)  call  SAVMATFORM 

stop 

end 


cl  SIMULATE:  subroutine  to  obtain  a  time  history  for  the  a/c  model  I 

c - - 

subroutine  SIMULATE(yth,u,nsteps) 
include 'MFBIDS.INC 

dimension  yth(maxout,2*maxtsteps),  u(2’''maxtsteps) 

double  precision  y(maxout) 

double  precision  atol,  rtol,  xstate(maxstates) 

double  precision  i-work(22*l()*maxstates+(2*l+l)*maxstates) 

double  precision  tstart.tend,  u0,ul 

integer  neq,itol,itask,istate,iopt,lrw,liw,mf 

integer  iwork(20+maxstates) 

common  /eqsmotcom/y,  tstart,  tend,uO,u  1 

c  **for  guidance  in  selecting  sizes  for  arrays  rwork  and  iwork 
c  see  source  code  LSODE.F 

external  EQSMOT 

c  **Setup  input  for  LSODE 

c  see  LSODE.F  for  explanation  of  these  parameters 

c 

neq=nstates 
c  tstart=starting  time 
c  tend=ending  time 
itol=l 

rtol=1.0E-06 

atol=1.0E-10 

itask=l 

istate=l 

iopt=0 

c  nvork=real  work  space 

lrw=22*  10*maxstates+(2*  1 + 1  )*max.states 
c  iwork=imiginary  work  space 
liw=20+maxstates 
c  iaceqs 
mf=23 

do  1  i=l,maxstates 
1  xstate(i)=0.() 


n  n 


do  11  i=l,maxout 
11  yth(i,l)=0.() 

tstart=0.0 

tend=0 

do  2  i=2,nsteps 
tend=tend+deltat 
uO=u(i-l) 
ul=u(i) 

call  LSODE(EQS  MOT.neq  ,xstate,tstart.tend,itol  ,rtol  ,atol , 
+  itask,istaie,iopt,rwoik,lrw,iwork,iiw,jaceqs,mO 
c  **If  istate  is  not  equal  to  2  then  an  error  has  occurred 
c  and  the  following  will  be  printed 

if  (istate.ne.2)  write(6,*)  'istate  = istate, 

+  'tend  = tend 

tstart=tend 
do  3  j=l,nout 
3  ythO,i)=y(i) 

2  continue 

return 

end 


I  DAT AIN;  subroutine  to  read  data  from  file  MFBIDS.INP  I 

c . 

c 

subroutine  DATAIN 
include 'MFBIDS.INC 
character*  80  ctitle 
character*  80  cease 
character*  15  edata,  cmatrixx 
common  /case/ccase 
common  /title/ctitle 
common  /fnames/edata,  cmatrixx 
character  cdummy*80 

c  **read  model  title 

read(unit=5,fmt='(a80)')  cdummy 
read(unit=5,fmt='(a80)')  ctitle 

c  **read  case  title 

read(unit=5,fmt='(a80)')  cdummy 
read(unit=5,fmt='(a80y)  cease 

c  **read  nstates,  nout,  noutmx 

read(unit=5,fmt='(a80)’)  cdummy 
read(unit=5,fmt=*)  nstates.nout,  noutmx 

c  **read  tmaximp,  deltat,  nsubintvl 
read(unit=5,fmt='(a80)')  cdummy 
read(unit=5,fmt=*)  tmaximp,  deltat 
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c  **read  sigma,  akmin,  akmax,  nkvals 
read(unit=5,fmt='(a8())')  cdummy 
read(unit=5,fmt=*)  sigma,  akmin,  akmax,  nkvals 

c  **read  iouttype,  impres,  ixcres,  iwave 
read(unit=5,fmt=’(a8())')  cdummy 
read(unit=5,fmt=*)  iouttype,  impres,  iexcres,  iwave 

c  **read  data  file  name 

read(unit=5,fmt='(a8())')  cdata 
read(unit=5,fmt='(al2)')  cdata 

4 

c  **read  matrixx  file  name 

read(unit=5,fmt='(a8())')  cdummy 
read(unit=5,fmt='(al2)')  cmatrixx 

close(unit=5) 

return 

end 

c . 

cl  SAVMATFORM:  subroutine  to  write  output  to  a  MATRIXx  readable  rile  I 

c . 

subroutine  SAVMATFORM 
include  'MFBIDS.INC 
dimension  rtemp(2*maxtsteps,maxstates) 
character*  10  vamame,  cnout 
character*  15  cdata,  cmatrixx 
common  /fnames/cdata,  cmatrixx 

open(unit=  1  ,file=cmatrixx) 
rewind(l) 

c  ** write  gust  intensity 
rtemp(l,l)=sigma 

call  MATSAV(l,'sigmag',  2*maxtsteps,  1, 

,  +  1, 0,  rtemp,  rtemp, '(Ip2e24.15y) 

c  **write  k  values  for  each  case 
»  do  1  i=l, nkvals 

1  rtemp(i,l)  =  allkvals(i) 

call  MATSAV(l,'kvals',  2*maxtsteps,  nkvals, 

+  1, 0,  rtemp,  rtemp,  '(Ip2e24.I5)') 

c  **write  which  output  quantity  was  "maximized" 
rtemp(l,l)=noutmx 

call  MATSAV(l,'noutmx',  2*maxtsteps,  1, 

+  1,0,  rtemp,  rtemp, '( 1  p2e24. 1 5)’) 

c  ** write  time  step 
rtemp(l,l)=deltat 

call  MATSAV(l,'deltat',  2*maxtsteps,  1, 

+  1,0,  rtemp,  rtemp, '( 1  p2e24. 1 5)') 
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c  **write  tmaximp 
rtenip(  1 .  l)=tmaximp 

call  MATSAV(1, 'tmaximp',  2*maxtsteps,  1. 

+  1,0,  rtemp,  rtemp,  '(Ip2e24.15)’) 

c  * '•‘write  maximum  output  values 
do  3  i=l,nout 
do  3  k=l,nkvals 

3  rtemp(k,i)=aexc_res(ntsteps,i,k) 

call  MATSAV(l,'maxout',  2*maxtsteps,  nkvals, 

+  nout,  0,  rtemp,  rtemp,  '(Ip2e24.15)’) 

c  **write  impulse  responses 
if  (impres.eq.l)  then 
do  4  i=l,nout 

open(unit=50,  status='SCRATCH') 
write(unit=50,fmt=  '(i3)')  i 
rewind(50) 
ndigits=l 

if  (i.gt.9)  ndigits=2 
if  (i.gt.99)  ndigits=3 
cnout='  ' 

read(unit=5{),fmt='(a3)')  cnout(  hndigits) 
close(50) 

vamame='impres'//cnout 
do  5  k=l, nkvals 
do  5  nt=l,ntsteps 

5  rtemp(nt,k)=aimp_res(nt,i,k) 

call  MATSAV(1,  vamame,  2*maxtsteps,  ntsteps, 
+  nkvals,  0,  rtemp,  rtemp, '( Ip2e24. 1 5)') 

4  continue 
endif 

c  ♦♦write  excitation  responses 
if  (iexcres.eq.l)  then 
do  6  i=i,nout 

open(unit=50,  status='SCRATCH') 
write(unit=50,fmt=  '(i3)')  i 
rewind(50) 
ndigits=:l 

if  (i.gt.9)  ndigits=2 
if  (i.gt.99)  ndigits=3 
cnout='  ' 

read(unit=50,fmt='(a3)')  cnout(  1  mdigits) 
close(50) 

varname='exresp'//cnout 

do  7  k=l, nkvals 
do  7  nt=l,(2^ntsteps-l) 

7  rtemp(nt,k)=aexc_res(nt,i,k) 
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call  MATSAV(1,  vamarne,  2*maxtsteps,  (2*ntsteps-l), 
+  nkvals,  0,  rtemp,  rtemp, '( Ip2e24. ISy) 

6  continue 
endif 

c  **write  excitation  waveforms 
if  (iwave.eq.l)  then 

open(unit=50,  status='SCRATCH') 
write(unit=5(),fmt=  '(i3)')  noutmx 
rewind(50) 
ndigits=l 

if  (noutmx.gt.9)  ndigits=2 
if  (noutmx.gt.99)  nciigiLs=3 
cnout='  ' 

read(unit=50,fmt='u3)')  cnout(l:ndigits) 
close(50) 

vamame='wavef//cnout 
do  9  k=l, nkvals 
do  9  nt=l,ntsteps 
9  rtemp(nt,k)=wave(nt,k) 

call  MATSAV(1,  vamarne,  2*maxtsteps,  ntsteps, 

+  nkvals,  0,  rtemp,  rtemp, '( Ip2e24.15)') 

endif 

close(l) 

return 

end 


I  SAVASCIIFORM:  subroutine  to  write  output  to  an  scii  file  I 

c . . . - 

subroutine  SAVASCIIFORM 
include  'MFBIDS.INC 
character*  80  ctitle 
character*  80  cease 
character*  15  edata,  cmatrixx 
common  /case/ccase 
common  /title/ctitle 
common  /fnames/edata,  cmatrixx 

open(unit=9,  file=cdata) 
rewind(9) 

c  **write  general  infomiation 
write(9,'(a80)')  ctitle 
write(9,'(a80y)  cease 
write(9,*) '  sigma  = ',  sigma 
write(9,*) '  matched  output  quantity  =  noutmx 
write(9,*) '  tmaximp  =  ',  tmaximp 
write(9,*) '  deltat  = ',  deltat 


! 


c  **write  maximum  output  values  for  each  case 
write(9.*) ' ' 

write(9  *)  '***************************************************' 
write(9’.*)  ‘MAXIMIZED  AND  TIME  CORRELATED  MAX  OUTPUT 
QUANTITIES' 

write(9,''‘)  '***************************************************' 
do  I  i=l,nout 

write(9,*) '  output  quantity  =  i 
write(9,*)  ’  k  value  maximum  output  value' 
do  I  k=I,nkvals 
write(9,*)  allkvals(k), ' 

+  aexc_res(ntsteps,i,k) 

1  continue 

c  **write  impulse  responses 

if  (impres.eq.  1)  then 
do  2  k=I,nkvals 
write(9,*) ' ' 

write(9,*)  '**+*************’•'***♦*****+***************♦****' 
write(9!*)  'IMPULSE  RESPONSES' 

write(9  *)  '********************************’*‘**************' 
write(9,*) '  kvalue  = ',  allkvals(k), 

+  '  maximized  output  quantity  = ',  noutmx 

do  2  iout=I,nout 

write(9,*)  'output  quantity  =  ',  iout 
do  2  j=I,ntsteps 
write(9,*)  aimp_res(j,iout,k) 

2  continue 
endif 

c  **write  excitation  waveforms 
if  (iwave.eq.l)  then 
do  3  k=l,nkvals 
write(9,*) ' ' 

write(9  *)  '*’*‘***+************************'*‘********’t‘*******' 
write(9’.*)  'EXCITATION  WAVEFORMS’ 

write(9,*)  '************+ ****************+*********** ******' 
write(9,*) '  kvalue  =  ',  allkvals(k), 

+  '  output  quantity  = noutmx 

do  3  j=l,ntsteps 
write(9,*)  wave(j,k) 

3  continue 
endif 

c  **write  excitation  responses 
if  (iexcres.eq.l)  then 
do  4  k=l,nkvals 
write(9,*) ' ' 

write(9  *)  '*’t‘***’t‘****’i‘****’t‘'t‘++******'*‘*******’i‘*’'‘***********' 
write(9!*)  'EXCITATION  WAVEFORM  RESPONSES’ 
write(9  *)  '**************’t"f‘*’*‘*******‘t'******’'‘*****’t‘******’i‘*' 
write(9,*) '  kvalue  =  allkvals(k). 
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+  '  maximized  output  quantity  =  noutmx 

do  4  iout=l,nout 

write(9,*)  'output  quantity  =  iout 
do  4  j=i,2*ntsteps-l 
write(9,*)  aexc_res(i,iout,k) 

4  continue 
endif 


close(9) 

return 

end 


MATSAV:  write  variables  to  a  file  in  matrixx  format  I 

c . . 

subroutine  MATSAV  ( lunit,  name,  nr,  m,  n,  img, 

+  xreal,  ximag,  formt ) 

c 

c  . . 

c 

c  MATSAV  writes  a  matrix  to  a  file  in  a  format  suitable  for  the 
c  matrixx  load  operation, 
c 

c  - 

c  parara.  type  on  input-  on  output- 

c  - 

c 

c  lunit  integer  fortran  logical  unit  number.  unchanged, 
c 

c  name  character*(*)  name  of  the  matrix,  one  al-  unchanged. 


c 

c 

(maximum 
length  10) 

phabetic  followed  by  up  to  9 
alphanumeric  characters. 

w 

c  nr 

c 

c 

c 

c 

integer 

row-dimension  in  the 
defining  dimension  or  type 
statement  in  the  calling 
program,  nr  must  be  greater 
than  or  equal  to  m. 

unchanged. 

c  m 

n 

integer 

number  of  rows  of  the  matrix 

unchanged. 

C 

c  n 

c 

r* 

integer 

number  of  columns  of  the 
matiix. 

unchanged. 

L 

c  img 
c 

c 

integer 

if  img  =  0,  the  imaginary 
part  (ximag)  is  assumed  to 
be  zero  and  is  not  saved. 

unchanged. 

c 

c  xreal 
c 

double 

precision 

real  part  of  the  matrix  to 
be  saved. 

unchanged. 

c  ximag  double 

imaginaiy  part  of  the  matrix 

unchanged. 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 


precision  to  be  saved. 

formt  characier*(*)  string  containing  the  for-  unchanged, 
(maximum  tran  format  to  be  used  for 
length  20)  writing  the  elements  of  the 
matrix. 


example:  the  following  fortran  program  generates  an  elementary 
matrix  in  x  and  writes  it  to  fortran  unit  1 .  assume 
that  unit  1  has  been  preallocated  as  file  (data  set) 
test. 


dimension  x(20,3),  dummy 
do  2(X)  j=l,3 
do  10()i=l,10 
x(i.j)=0.0d0 
100  continue 
x(j,j)=1.0d0 
200  continue 

call  MATSAV(  1,  'amatrix',  20,  10,  3, 0, 
$  X,  dummy,  '(Ip2e24.15)’ ) 

stop 
end 


after  this  program  runs,  invoke  matrixx  and  type: 

<>  load  'test' 

this  will  put  x  on  the  stack  as  stack- variable-name  amatrix. 


integer  lunit,  m,  n,  nr,  img 
character*(*)  name,  fount 
dimension  xreal(nr,l),  ximag(nr,l) 
character  nam*10,  form*20 


write  header  record. 


nam=name 

fonn=formt 

write(lunit,'(a  1 0,3i5,a20)')  nam ,m ,n,img,form 


write  real-part  of  the  matrix. 


write(lunit,form)  ((xreal(i,j),i=l,m),j=l,n) 


write  imaginary-part  if  nonzero. 
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if(img.ne.O)  write(lunit,fonn)  ((ximag(i,i),i=l,ni),j=l,n) 

return 

end 


« 


t 
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Appendix  B  •  MFBIDS.INC  source  code  and  description  of  parameters 

parameter(maxtsteps=5(X)l,  maxstates=4(),  maxout=20,  maxkvals=20) 
common  /modelinfo/nstates,  nout 
common  /siminfo/tmaximp,  deltat,  ntsteps,  nsubintvl 
common  /case_info/noulmx 

common  /excitation_intb/sigma,  akmin,  akmax,  nkvals, 

+  allkvals(maxkvals) 

common  /t_hist/aimp_res(maxtsteps,maxout,maxkvals), 

+  aexc_i'es(2*maxtsteps,maxout.maxkvals), 

+  wave(maxtsteps.maxkvals), 

+  aimpt(2*maxt.steps) 

common  /output/iouttype,  impres,  iexcres,  iwave 


The  following  is  a  description  of  the  parameters: 

maxtsteps  This  parameter  is  the  maximum  number  of  time  steps  for  the 

impulse  responses,  (maxtsteps  >  tmaximp/deltat  +1) 

maxstates  This  parameter  is  the  maximum  number  of  states  in  the  equations 

of  motion,  (maxstates  >  nstates) 

maxout  This  the  maximum  number  of  output  quantities  required  by  the 

aircraft  model,  (maxout  >  nout) 

maxkvals  This  the  maximum  number  of  kvalues  which  can  be  run. 

To  change  these  parameters  edit  MFB  IDS.INC.  In  addition,  maxout  is  also  found  in 
MODEL.F  and  must  be  equal  that  the  value  in  MFBIDS.INC.  Computer  storage  can  be 
minimized  by  making  maxstates  and  maxout  equal  to  the  minimum  values  required  for  a 
given  analytical  model. 


i 


f 
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Appendix  C  -  MODEL.F  (ARW-2)  source  code  listing 


subroutine  EQSMOT(neq,t,x,xd) 
parameter(maxout=2()) 
double  precision  x(*),  xd(*),  y(maxout) 
double  precision  t.tstart.tend,  uO,  u  1 
common  /eqsmotcom/y  ,Lstart.tend.u(),u  1 

u=(t-  tstart)  *  (u  1  -  u())/(tend-  tstart)+uO 


c . - . state-space  system 

i  c  —  {nlarw2. controller.  1 } 

C  ““  SS  C 

yd)  =  -5.24497()395787696D-4*x(l)  -  2.747(X)46914()994399D-4*x(2) 

•  +  -  1.42946149773359D-3*x(3)  +  5.56258182659891597D-6*x(4) 

y(2)  =  -3.287044785655 141()2D-3*x(l)  -  6.743878267402364D-5*x(2) 
+  -  9.87258 1358599253()6D-4*x(3)  + 

-i-  9.90556693932225  l()4D-6*x(4) 

c . - 1  -  u  bounded  limit 

c  —  {nlarw2.aileron  limiter.2} 

y(3)  =  min(  O.ODO,  max(  -0.0174499999999999999D0,  y(2) ) ) 
c  y(3)  =  min(  0.0 1 74499999999999999D0,  max( 
c  +  -0.0174499999999999999D0,  y(2) ) ) 

c . state-space  system 

c  ”  {nlarw2.arw2  moD.3} 

Q  --  SS  c  •• 

y(4)  =  -2846.47376 1074609D0*x(5)  -  2699.5767200 1826D0*x(6)  + 

+  6978.2941 1601534196D0*x(7)  +  2. 14350293 155(X)8902D5*x(8)  - 

+  22.3289763900712601D0*x(9)  -  3.3990629437 14242D0*x(  10)  - 

+  9.13l44768834393994DO*x(l  1)  +  5.05173942454337099D0*x(12) 

+  +  75.7538  i09829206405D0*x(l3)  - 

+  0.95882 1750950846804D0*x(14)  - 

+  0.86024628749677 1303D0*x(  15)  -  78.95852907402 15004D0*x(l 6) 

-H  +  43.922 1734285092698D0*x(  17)  + 

+  8.35298937508213202D0*x(18)  -  0.95882 1750950846804D0*x(  19) 

+  -  0.86024628749677 1 303D0*x(20) 

y(4)  =  y(4)  -  78.95852907402 1 5(X)4D0*x(21)  + 

-t-  43.922 1734285092698DO*x(22)  +  8.352989375082 13202D0*x(23) 

^  +  -  1954.93624505O99399D0*x(24)  - 

+  2.4341 670262402 1802D0*x(25)  -  2597.51251 185656298D0*x(26) 

,  +  +  37578.7624927093302D0*x(27)  - 

+  4.17208401208873203D0*x(28)  +  3.077920932(X)968897D9*x(29) 

+  -  41.0633443973176799DO*x(35)  + 

+  3.573 16720948225297D-3*x(36) 

y(5)  =  -42.1666883238785903DO*x(5)  +  I6.36587051303695D0*x(6) 

+  -  453. 129586 109553202DO*x(7)  +  35 1.46043821 35 13(X)lDO*x(8) 

+  -  0.59975463706 1799697DO*x(9)  - 

+  4.63233083 144964902D-2*x(  10)  -  1.5201 16146415901D-2*x(ll) 

-  6.86370728184497701D-2*x(12)  + 

+  0.48369850553l560102D0*x(13)  - 

+  0.969166155996266099DO*x(14)  + 

+  0.519949869945893497D0*x(15)+  l.26302828846864701D0*x(16) 

+  -  1. 7400070959 1828(X)lDO*x(  17)- 

+  1.1 89239579268 164D-2*x(  18)  -  ().969166155996266099D0*x(19) 
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+  +0.519949869945893497DO*x(20) 

y(5)  =  y(5)  +  1. 26302828846864701  D0*x(21)  - 
+  1.7400070959i828001D0*x{22)-  1.189239579268  164D-2*x{23) 

+  +  900.378460475793005DO*x(24)  + 

+  0.64379949(X)67852503D0*x(25)  +  1238.1447(XX)8882901DO*x(26) 

+  +  291. 10767237057 1898DO*x(27)  - 

+  9.35961201 178221493D-2*x(28)  +  3.81608823019032497D7*x(29) 

+  -  0.6521 16307 184986296DO*x(35)  - 

+  2.32409992594539401D-5*x(3o) 

y(6)  =  6384.77 199999999698D0*x(6)  -  2628 1.1 30000000000 lD0*x(7) 

+  -  44252.8600(XXXXX)997D0*x(8) 

y(7)  =  81.7875899999999092D0*x(6)  -  1061.036D0*x(7)  - 
+  2228.2‘^D0*x(8) 

y(8)  =  -231.74300000(XXX)4D0*x(6)  -  3091.725999999995D0*x(7)  + 

+  61 62.58 10{XXXXXX)604D0*x(8) 

y(9)  =  28.7296599y99999702DO*x(6)  +  534.1 0220000000 1204D0*x(7) 
+  -  579 1 .63599999999894D0*x(8) 

y(10)  =  x(24) 
y(ll)  =  x(25) 
y(12)  =  x(27) 
y(13) =  x(28) 
y(14)  =  x(35) 

y(15)  =  -2846.47376 1074609DO*x(5)  -  2699.5767200 1826D0*x(6)  -h 
+  6978.2941 1601534196D0*x(7)  -i-  2. 14350293 155008902D5*x(8)  - 

+  22.328976390071 260 lD0*x(9)  -  3.3990629437 14242D0*x(10)  - 

+  9.13144768834393994D0*x(l  1)  +  5.05173942454337099D0*x(12) 

-t-  +75.7538109829206405D0*x(13)- 

+  0.95882 1750950846804D0*x(14)  - 

+  0.86024628749677  L303D0*x(  15)  -  78.95852907402 15004D0*x(  16) 

+  43.922 1734285092698D0*x(  17)  + 

+  8.352989375082 13202D0*x(l  8)  -  O.958821750950846804D0*x(19) 

-I-  -  0.86024628749677 1 303D0*x(20) 

y(15)  =  y(15)  -  78.9585290740215004D0*x(21)  + 

-h  43.922 1734285092698D0*x(22)  +  8.352989375082 13202D0*x(23) 

+  -  1954.93624505099399D0*x(24)  - 

-I-  2.4341670262402 1802D0*x(25)  -  2597.51251 185656298D0*x(26) 

+  +  37578.7624927093302D0*x(27)  - 

4.17208401208873203D0*x(28)  +  3.07792093200968897D9*x(29) 
+  -  41.0633443973 176799D0*x(35)  + 

+  3.573 16720948225297D-3*x(36) 

y(16)  =  -42.1666883238785903DO*x(5)  +  16.3658705 1303695D0*x(6) 
+  -  453.129586109553202D0*x(7)  +  351.460438213513001D0*x(8) 

-t-  -  0.59975463706 1799697D0*x(9)  - 

+  4.63233083 144964902D-2*x(  10)  -  1.5201 16146415901D-2*x(ll) 

-1-  -  6.86370728 18449770  lD-2*x(  12)  + 

+  0.48369850553 1560102D0*x(13)  - 

+  0.9691661.‘=5996266099D0*x(14)  + 

+  0.519949869945893497D0*x(15)+  1.26302828846864701D0*x(16) 

+  -  1.7400070959 182800  lD0*x(  17)- 

+  1. 189239579268 164D-2*x(  18)  -  0.9691 66 155996266099D0*x(  19) 

+  +  0.519949869945893497D0*x(20) 

y(16)  =  y(16)  +  1.26.302828846864701D0*x(21)  - 
+  1.740(X)709591828001D0*x(22)-  1.189239579268164D-2*x(23) 

+  +  900.378460475793(X)5D0*x(24)  + 
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+  0.64379949{X)678525()3D()*x(25)+  l238.1447(XK)8882901D0*x(26) 

+  +  291. 10767237057 1898D()*x(27)  - 

+  9.35961201 17822 1493D-2*x(28)  +  3.81608823019032497D7*x(29) 

+  -  0.6521 16307184986296D0*x(35)  - 

+  2.32409992594539401  D-5*x(36) 

c . 1 '  u  bounded  limit 

c  "  {nlarw2.evel  limiier.4} 

y(17)  =  min(  0.0174499999999999999D0.  max( 

+  -0.0 1 74499999999999999D0,  y ( 1 ) ) ) 

c . state-space  system 

c  —  {nlarw2.controller.l } 


c  —  ss  ab  " 

xd(l)  =  -27.797 1207546241899D0*x(l)  -  25.4156325927014D0*x(2) 

+  +  0.96687429703 1 38 1 502D0*x(3)  - 

+  8.71552238769890408D-2*x(4)  -  9.51654087303377394D-3*y(4) 

-  1.6803237 14059 14D-2*y(5) 
xd(2)  =  -6.876026223639.34302D0*xri)  - 

-t-  7.3671 1056293708599D0*x(2)  -  8.27135760534558596D-3*x(3) 

-t-  0.976 139521 572 193303D0*x(4)  - 

+  2.63798849248624(X)lD-3*y(4)-i-  1.501039866773363D-2*y(5) 

xd(3)  =  -52.3804248827 116 199D0*x(l)  - 
+  35.58994104867247()lDO'^x(2)  -  0.29(X)86695207493903D0*x(3) 

-  0.218148213265887DO*x(4)  -  7.09094969390322604D-3*y(4) 

+  -  9.45879838024 1 606()3D-3*y(5) 

xd(4)  =  -4067.66699224439799D0*x(l)  - 
-t-  5071. 334 17368828702D0*x(2)  -  5.043376280250499D0*x(3)  - 

+  17.0253581968938801D0*x(4)-  1.15740514598378101D-2*y(4) 

+  +  1. 73620839739 127D-2*y(5) 

c . state-space  system 

c  "  {nlarw2.arw2  moD.3} 
c  "SS  ab" 

xd(5)  =  x(10) 
xd(6)  =  x(  1 1 ) 
xd(7)  =  x(12) 
xd(8)  =  x(13) 

xd(9)  =  -1 25.2459238677.3790 lD0*x{5)  - 

45.03350413 19353904D0*x(6)  -  1 17.872 130364775799D0*x(7)  -i- 
+  3906.6306 1647083703D0*x(8)  -  1. 26414087776252901  D0*x(9)  - 

+  0.1273297457868026D0*x(10)  -  0.2001582942733044D0*x(l  1)  -h 

+  0.1463652541766454D0*x(12^+  l.l5969827063154399D0*x(13) 

-t-  -  0.969152580431533295D0*x(14)  -»■ 

-I-  8.27755253907873397D-4*x(15)  + 

+  6.04388965594990202D-3*x(16)- 

+  6.85614516592095201D-3*x(17)  + 

-t-  5.01807022310764295D-4*x(18)  - 

+  0.969 1 5258043 1533295DO*x(l 9) -t- 

8.27755253907873397D-4*x(20) 
xd(9)  =  xd(9)  6.04388965594990202D-3*x(21)  - 
+  6.85614516592095201D-3*x(22)  + 

+  5.01 8070223 10764295D-4*x(23l  +  1223.06162596226201D0*x(24) 

-t-  0.899244l09822731702D0*x(25)  + 

1763.75969 19 1059301D0*x(26)  +  698. 16 1055555356  lDO’^x(27)  - 
-t-  0. 1264279353258893DO*x(28)  -(■  6.7 1666950223660506D7*x(29) 

-»■  -  1. 867833656 108772D0*x(35)  +  7.6556232257 1364598D-5*x(36) 


c-3 


+  -  4.845331 15551499498D-4+U 

xd(lO)  =  -3 1.7973764346982  l{)2D()*x(5)  -  16.602 17980169()8D0*x(6) 

+  -  135.2985761432064D0*x(7)  +  2147.94963514691301D0*x(8)  - 

+  1.24869739759363()99D0*x(9)  -  0.2655459480335534D0*x(10)  - 

+  0.153253792353896599DO*x(l  1)  -  1.0508588744 17987D-2*x(  12) 

+  +  1.278272593746799D0*x(13)  -  1. 27530495555 126799D-3*x(  14) 

+  -  5.79125135412297698D0*x(15)  + 

+  2.918827546093392D-3*x(16)-  1.20106974790340501D-3*x(17) 

+  +  3.6673 133776078 12D-3*x(  18)- 

+  1.27530495555126799D-3*x(19)  -  5.79 1251 354 12297698D0*x(20) 

+  +2.918827546093392D-3*x(21) 

xd(lO)  =  xd(lO)  -  1.20106974790340501D-3*x(22)  + 

+  3.6673 133776078 12D-3*x(23)  +  4267.34662588543097D0*x(24) 

+  +  3.5 1 37358075365630 1  D0*x(25)  + 

+  6539.00920066869(K)3D0*x(26)  +  779.84733353379 12D0*x(27)  - 

+  0. 10065335997 10573D0*x(28)  +  5.69431331591210403D7*x(29) 

+  -  0.53085731 374576990  lD0*x(35)  + 

+  6.99578 120780643495D-5*x(36)  -  4.427709625 193978D-4*u 

xd(l  1)  =  -4035.453860038237DO*x(5)  -  5039.804988 19288903D0*x(6) 
+  -  103 15.3494347062699DO*x(7)  +  2.4223448237 1904D5*x(8)  - 

+  27.9029350865074499D0*x(9)  -  5.(K)483961999364602D0*x(10) 

+  -  16.9199342752880302D0*x{l  1)  - 

+  1.21 7457408  l8824499D0*x(  12)  +  75.01 44208438646292D0*x(  13) 

+  +  8.243(K)236275332708D-3*x(  14)  + 

+  1.93965721381322799D-3*x(15)-  152.1 134291940743D0*x(16) 

+  -  5.8038654 1425226303D-2*x(  17)  + 

+  5.25523764977082796D-2*x(18)  + 

+  8.24300236275332708D-3*x(19)  + 

+  1 .9396572 1381 322799D-3*x(20) 

xd(  1 1 )  =  xd(  1 1 )  -  1 52. 1 1 3429 1 940743D0*x(2 1 )  - 
+  5.80386541425226303D-2*x(22)  + 

+  5.25523764977082796D-2*x(23)  -  3324.495306375349D0*x(24) 

-H  -  3.661 54656 177599203D0*x(25)  - 

+  3889.2879735487080  lD0*x(26)  +  70121.2766765817105D0*x(27) 

+  -  16.230958826885O602DO*x(28)  + 

+  8.2431213221 121 5 198D9*x(29)  -  58.8 1479488 18334498D0*x(35) 

+  +  6.03774882724283902D-3*x(36)  - 

+  3.8213600172423 1499D-2*u 

xd(12)  =  387 1.7533554938 1 299D0*x(5)  +  1535.837069757428D0*x(6) 
+  -  34825.31 67988588 102D0*x(7)  -  1.03428353574103699D5*x(8) 

+  +  26.295468997 1899098D0*x(9)  +  1.250989170291227D0*x(10) 

+  +  2.630979503694448D0*x(l  1)  -  18.3240096889994701D0*x(12) 

+  +  9.72678577610167805D0*x(13)  - 

+  4.201 1207881 598 1804D-3*x(  14)  + 

+  1.3530565 1065 17070  lD-3*x(  15)  + 

+  4.063760 124791 73 199D-3*x(  16)  -  140.368755 107959299D0*x(  17) 

+  +  0.141684134413155902DO*x(18)  - 

+  4.201 12078815981804D-3*x(19)  + 

+  1 .35305651065 1 7070 1  D-3*x(20) 

xd(12)  =  xd(12)  +  4.063760 12479 173 199D-3*x(21)  - 
+  140.368755 107959299D0*x(22)  +  0.1416841344 13 155902D0*x(23) 

+  2449.12248315986699D0*x(24)  + 

+  2.20529981 83453 1003D0*x(25)  +  1948.391 10573395001D0*x(26) 

+  -H  20503.2354698 1 1 2499D0*x(27)  - 
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+  8.717721341283493()8D()*x(28)  +  3.55976()45603601098D9*x(29) 

+  +  56.3946 108862526 199D()*x(35)  - 

+  3.633780767338363()2D-3*x(36)  +  2.29986124515086901D-2*u 

xd(13)  =  -463.806180072810l01DO*x(5)  - 
+  349.736890131978697D0*x(6)  -  1 195.977582526502D0*x(7)  - 

+  40097.8732828805196D()*x(8)-  1.46494257 127 194799D0*x(9)  + 

+  0.1 83876792527 185799D0*x(  10)  - 

+  0.5551 17409039766599D0*x(l  1)  + 

+  0.167101544533691599D0*x(12)  -  29.83280229208 18499D0*x(13) 

+  -  3.81290566534986402D-3*x(14)  + 

+  5.78150706319158894D-4*x(15)  -  0.1002548480818097D0*x(16) 

,  +  +  9.99420602 1040908D-2*x(17)  -  6.49269676792 184203D0*x(l  8) 

+  -  3.81290566534986402D-3*x(19)  + 

+  5.78150706319158894D-4*x(20) 

»  xd(13)  =  xd(13)  -  0.1002548480818097D0*x(21)  + 

+  9.99420602 1040908D-2*x(22)  -  6.49269676792 184203D0*x(23) 

+  +  1019.796982 139836D0*x(24)  +  0.9769050070591 09 196D0*x(25) 

+  +  2093.9417 1462465602D()*x(26)  - 

+  5232.10908(K)6390696D0*x(27)  -  1.29605906458 1286D0*x(28)  + 

+  1. 28570 1909877968D8*x(29)  -  7.1 1996318155942698D0*x(35)  + 

+  6. 163729505 18372805D-4*x(36)  -  3.90109462353402403D-3*u 

xd(14)  =  0. 12848 1999999999999D0*x(9)  - 
+  40.1736(X)0(KKM)00804DO*x(10)  +  4.33304999999998597D0*x(l  1) 

+  -  98.3560999999999694D0*x(12)  +  1044.980000000003D0*x(13) 

+  -  84.254265530464 1794D0*x(  14)  - 

+  146.22500(KKXKKK)399D0*x(25)  +  158.846999999999799D0*x(28) 

+  +  8.527042966739 1 8397D-2*x(35)  + 

+  4.06559376838870898D-2*x(36)  -  0.25731 606 12904246D0*u 

xd(15)  =  -19.28531KXKK)00{K)101D0*x(9)  - 
+  5.5699 1999999999599D0*x(  10)  +  10.0029999999999899D0*x(  1 1) 

+  -29.1312000000(KX)402DO*x(12)- 

+  92.8560999999999694D0*x(  1 3)  -  84.254265530464 1794D0*x(  15) 

+  -  1 1 1 .353(H)(MX)(MX)0099D0*x(25)  + 

+  17.2043999999999602D0*x(28)  +  4.65810662194168197D-2*x(35) 

+  +  2.2209304360874 1 999D-2*x(36)  - 

+  0.1405652 1747388740  lD0*u 

xd(16)  =  6.8972  l(K)0(KKKXKK)95DO*x(9)  - 
+  6.3203900(KXX)(XK)295D0*x(10)  -  5.16849999999999499D0*x(l  1) 

^  +  +  2.07699999999999801D0*x(12)  + 

+  205.770999999999699D0*x(13)  -  84.2542655304641794D0*x(16) 

+  +  1 9  065 1999999999999D0*x(25)  - 

*  +  59.3106999999999998D0*x(28)  +  3.427498257283779D-4*x(35) 

+  +  1.634190845563998D-4*x(36)  -  1.034298(X)3521519D-3*u 

xd(17)  =  -7.4448 19999999993DO*x(9)  + 

+  6.893470(XX)0(KXX)798D0*x(10)  +  8.567650000{XXX)1509D0*x(l  1) 

+  +  2.0852D0*x(12)  -  145.588999999999899D0*x(13)  - 

+  84.254265530464 1794D0*x(  17)  -  15.7273(XXXXXXXX)lD0*x(25)  - 

+  32.4503999999999504D0*x(28)  +  3. 19475760574534596D-2*x(35) 

+  +  1.52322284 1035782D-2*x(36)- 

+  9.64065089263 157499D-2*u 

xd(18)  =  - 13.323899999999980 lD0*x(9)  - 
+  18.1652(XXXXXXXX)302D0*x(I0)-  118.12700(XXXXXXXX)lD0*x(ll) 

+  +  1 34.6 1 4999999999799D0*x(  1 2)  - 

+  949.75400(XXXX)00801D0*x(13)  -  84.254265530464 1794D0*x(  18) 
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+  +  3. 1790500()00000()398D()*x(25)  - 

+  397.355999999999803D0*x(28)  +  0.27 8893228004259  lDO*x(35) 

+  +0.132973010015630801DO*x(36)- 

+  0.841601329212856797D0*u 

xd(19)  =  82.22080000(XXXX)503D0*x(9)  + 

+  38.1446000000000804D0*x(10)  -  9.01630999999997607D0*x(l  1) 

+  +31.3315999999999799D0*x(12)  + 

+  656.121999999999403D0*x(13)  -  168.6879863655722DO*x(19)  + 

+  431.76300(XXXXXX)80IDO*x(25)+  123.543000000000101D0*x(28) 

+  +  0.7835 19507974698798D0*x(35)  + 

+  0.373572883525770998D0*x(36)  -  2.3643853387707 1699D0*u 

xd(20)  =  2.010990000(XXXX)699D0*x(9)  + 

+  29.2762000(XXXXX)2()2DO*x(10)  -  12.37 16999999999799D0*x(l  1) 

-«■  +  20.6739(XXXXXXXXXX)2D0*x(  1 2)  + 

+  247.386999999999698D0*x(13)  -  l68.6879863655722DO*x(20)  + 

+  264.54900(XXXXXX)902D0*x(25)  +  26.5588999999999902D0*x(28) 

+  +  1. 0063890730269250  lD-2*x(35)  + 

+  4.79834470147791304D-3*x(36)  -  3.03692702625 1846D-2*u 

xd(21)  =  22.594(XXXXXXXXX)5D0*x(9)  -  5.37708000000000697D0*x(10) 
+  +  8.3899499999999989  lD0*x(l  1)  -  18.90300000000002D0*x(12) 

+  +  142.493(XXXXXXXX)4D0*x(I3)  -  168.6879863655722D0*x(21) 

+  -  22. 1969(XXXXXXXX)30 1  D()*x(25)  + 

+  168.935999999999702D0*x(28)  +  0.223846017964695498D0*x(35) 

+  +  0. 106727 147883(X)2501D0*x(36)  - 

+  0.675488277740527096D0*u 

xd(22)  =  -23.1385000(XXXXX)201D0*x(9)  + 

+  4.77663999999998601D0*x(10)  -  12.28980000000001D0*x(l  1)  + 

+  13.8938999999999699D0*x(  12)  -  266.941999999999098D0*x(13) 

+  -  168.6879863655722DO*x(22)  +  19.15959999999995D0*x(25)  - 

+  43.0659(XXXXXXXX)596D()*x(28)  -  0.231364769970956501D0*x(35) 

+  -0.1 103 12(XX)38368 19()lD()*x(36)  + 

+  0.698177217618241(X)3D()*u 

xd(23)  =  150.079999999999899D0*x(9)  - 
+  14.4340999999999999D()*x(10)  +  184.738999999999599D0*x(l  1) 

+  -  262.1530(XXXXXXX)198D{)*x(12)  + 

-i-  2774.57000(XKXKX)698D()*x(13)  -  168.6879863655722D0*x(23)  + 

+  28.3 1230(XXXXXXX)5O2D0*x(25)  +  935.655999999998997D0*x(28) 

+  +  0.302353929991241402D0*x(35)  + 

+  0.144158796714765301D()*x(36)  -  0.91 239744756 1808406D0*u 

xd(24)  =  x(25) 

xd(25)  =  -3.94784 179999999702D5*x(24)  - 
-I-  1256.637000(XXXX)199D0*x(25)  +  7.89568359999999404D6*x(26) 

xd(26)  =  -20.0D0*x(26)  -t-  y(17) 
xd(27)  =  x(28) 

xd(28)  =  -3.348D5*x(27)  -  818.4(KXXXXXXXX)1498D0*x(28)  + 

+  7.201548D+ll*x(29) 

xd(29)  =  x(30) 

xd(30)  =  -2.151D6*x(29)  -  54().()99999999998502D0*x(30) 

+  2.191D5*x(31) 

xd(31)  =  x(32) 

xd(32)  =  -2.191D5*x(31)  -  185.3(XXXXXXXX)002D0*x(32)  + 

+  3.742D6*x(33) 

xd(33)  =  x(34) 

xd(34)  =  -3.742D6*x(33)  -  1446.5D0*x(34)  +  y(3) 
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xd(35)  =  -0.33999999y9999999D()*x(35)  - 
+  0.1621079999999999()lD()*x(36)+  1.0260(XXXX)000(X)299D0*u 

xd(36)  =  -0.359999999999999397D0*x(36)  +  u 

c - 

c 

return 

end 


t 


« 
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Appendix  D  •  MODEL.INP  file  for  the  ARW-2  model 


This  is  an  example  of  the  input  file.  MODEL.INP,  that  must  be  created  by  the  user.  The 
odd  numbered  lines  are  dummy  variables  used  to  name  the  fields  and  the  even  numbered 
lines  contain  the  actual  data.  This  is  the  input  file  used  to  generate  the  results  shown  in 
figure  5. 

model 

arw2  flexible  nonlinear  model,  much  .86,  altitude  24k  ft 
case 

maximize  output  #6 
nstates  noutnoutmx 
36  17  6 

tmaximp  deltat 
10.0  0.005 

sigma  akmin  akmax  nkvals 
1530.0  10.0  15000.0  9 
iouttype  impres  iexcres  iwave 
3  111 

standard  filename  (iouttype=l) 
arwl530.out 

matrixx  filename  (iouttype=2) 
arwl530.mat 


Where  these  input  quantities  are  defined: 

model  A  character  variable  describing  the  model  used  in  the  analysis 

case  A  character  variable  describing  the  case 

nstates  Number  of  states  in  the  model 

nout  Number  of  output  quantities  used  in  the  model 

tmaximp  The  length  of  the  simulation  used  in  the  impulse  responses 

deltat  The  time  step  used  in  the  simulations 

noutmx  The  number  of  the  output  quantity  to  be  maximized 
sigma  The  gust  intensity  used  in  the  analysis,  units  of  velocity 

akmin  The  minimum  k  value  to  be  used  in  the  analysis 

akmax  The  maximum  k  value  to  be  used  in  the  analysis 

)  nkvals  The  total  number  of  k  values  to  be  used  in  the  analysis 

1,  akmin  is  used 

2,  akmin  and  akmax  are  used 

f  >3,  k  values  distributed  logrithimically  between  akmin  and  akmax 

iouttype  Specifies  the  format  for  the  output 

1,  Normal  ascii  file 

2,  MATRIXx  readable  ascii  file 

3,  Both  1  and  2 

impres  If  1,  write  the  impulse  responses  to  the  output  file(s) 

iexcres  If  1,  write  the  excitation  waveform  responses  to  the  output  file(s) 

iwave  If  1,  write  the  excitation  waveform(s)  to  the  output  file(s) 
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Appendix  R  -  A  sample  listing  of  the  MFRIDS  program  output  using  the  input  file 
shown  in  apppendix  I)  and  the  model  listed  in  Appendix  C 

iir\v2  ricxibic  nonlincai  mocit'l.  Mach  86.  Altitude  24k  ft 
inaxiinizo  (uiff)iii  #6 

length  of  siimilaliun  for  impulse  responses  =  lO.OtKKH) 
mimhcr  of  time  steps  for  inpulse  responses  =  2(K)I 

length  of  simulation  for  cxcitalicm  responses  =  2().(KMK) 

numher  of  time  steps  for  excitation  responses  =  4()0I 

gust  intensity  =  l.'^.^O.OO 

calculating  impulse  responses  for  each  k  value 
calculating  noi  juali/ed  excitation  waveforms 


1 

k  = 

lO.OtKMM)  sqrtt energy )  (rf  output  6  = 

568.177 

T 

k  = 

24.9466  sqrtlcnergy)  of  output  6  = 

1417.29 

3 

k  = 

62.2333 

s(]rt(cncrgy)  of  output  6  = 

3536.37 

4 

k  = 

l.‘!.‘^.2.‘>l 

S(|rl(encrgy)  of  output  6  = 

8820.35 

5 

k  = 

387.298 

sc]rt(encrgy )  of  output  6  = 

22003.6 

6 

k  = 

966. 1 77 

s(|rl(cncrgy)  of  output  6  = 

56134.6 

7 

k  = 

2110.28 

si|rt( energy)  of  output  6  = 

162952. 

8 

k  = 

6012.84 

sqrltenci  gy)  of  output  6  = 

509979. 

9 

k  = 

1. ''000.0 

sc|rt(encrgy  )  of  output  6  = 

1.4941  ll-:+06 

calculating 

maximimi  responses 

1 

k  = 

lO.OOOOO 

maximum  value  of  output 

6  =  287()(H). 

2 

k  = 

24.9466 

tuaximum  value  of  output 

6  =  286965. 

3 

k  = 

62,2333 

maximum  value  of  output 

6  =  286988. 

4 

k  = 

15.5,251 

maximum  value  of  output 

6  =  286997. 

5 

k  = 

387.298 

maximum  value  of  output 

6=  287025. 

6 

k  = 

966. 1 77 

maximum  value  of  output 

6  =  289885. 

7 

k  = 

2410.28 

maximum  value  of  output 

6  =  296994. 

8 

k  = 

6012.84 

maximum  value  of  output 

6  =  279944. 

<) 

k  = 

I5(K)0.0 

maximum  value  of  output 

6  =  249730. 

saving  data 

•O.S.  GOVERNMENT  PRINTING  OPPICE:  199‘l-50'‘-078-00004 


